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CHAPTER 1 


INTRODUCTION 

Mathematics is a universal tool in the physical sciences, and 
much of the insight gained in other fields, especially in systems 
engineering, is directly applicable to hydrology (Dawdy, 1969). 

Since Modern Control and Estimation Theory have been applied success- 
fully to aerospace engineering problems (i.e., satellite tracking, 
orbit determination, space navigation, etc...) in the last two 
decades, and there exist many similarities between satellite 
tracking problems and the identification of unknown parameters of 
hydrologic processes (i.e., the models are not known precisely; the 
system under study is stochastic and highly non-linear; also, there 
is noise in the observations) , an attempt to use this new approach 
for the study of hydrologic systems is worthwhile to be investigated. 

Characteristics of Hydrologic Systems 

A hydrologic system may be defined as an interconnection of 
physical elements which are related to water in its natural state. 

The essential feature of a hydrologic system lies in its role in 
generating outputs (i.e., runoff,...) from inputs (i.e., rainfall, 
snowmelt, temperature,...), or in interrelating inputs and outputs. 
The stochastic nature of the inputs and outputs of hydrologic 
systems has been discussed by Yevjevich (1971) . 

Hydrologic processes are complex time-varying distributed 
phenomena, which are controlled by an unknown number of climatic 
and physiographic factors. The later descriptors tend to be static 
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or change slowly in relation to the time scale of hydrologic 
fluctuations. Observations of results in the laboratory and in the 
field also indicate that many of the component processes in hydrol- 
ogy are nonlinear due to the following reasons (Amorocho and Orlob, 
1961 ) : 

(1) the time variability of watersheds due to the natural 
processes of weathering, erosion, climatic changes, etc... 

(2) the uncertainty with respect to the space and time 
distribution of the inputs and outputs of hydrologic systems, and 
with respect to the states and characteristics of the interior 
elements of the system in time; and, 

(3) the inherent nonlinearity of the processes of mass and 
energy transfer that constitute the hydrologic cycle. 

Thus, for systems engineers, hydrologic processes can be considered 
as nonlinear dynamic distributed-parameter systems with partially 
known or unknown structures, operated in a continuously changing 
environment. The inputs and outputs of these systems are measurable 
but the data obtained are imbedded in noises with partially known 
or unknown characteristics. 

Modeling Problems in Hydrology 

Precise mathematical models developed for the study of 
hydrologic systems should be nonlinear, dynamic, distributed 

parameter models. However, at the present time, because of lack of 
data on parameter distribution, the assumption of space invariance is 
unavoidable. The subdivision of large watersheds into environmental 
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zones, where environmental conditions which affect the behavior 
of hydrologic systems can be assumed as uniform, and the use of 
lumped-parameter model for each zone, are then required to improve 
the modeling situation. By routing the flow spacewise through all 
the lumped— parameter models representing the environmental zones 
the total simulation of the entire watershed represents a distributed- 
parameter system. 

Lumped- parameter models of hydrologic systems can be divided 
into deterministic and stochastic models. The deterministic approach 
often is called parametric modeling. The choice of the model is 
determined by the type of problem to be solved. Parametric models 
require input data with considerable detail in time, therefore, they 
model transient responses well and are most widely used for short- 
term simulation or for actual prediction for water management pur- 
poses (Dawdy, 1969). Stochastic models have the advantage of 
taking into account the chance dependent nature of hydrologic events. 
Stochastic synthesis models are concerned with the simulation of the 
relationship between input and output data (cross-correlation models) 
and between successive values of each data series (serial correla- 
tion models). In stochastic simulation models, statistical measures 
of hydrologic variables are used to generate future events to which 
probability levels are attached. But in this case long term records, 
which in many instances are not available, are needed to estimate 
the parameters of the stochastic model in order to obtain a proper 
representation of their stochastic nature. Stochastic simulation 
models usually are used for planning purposes to develop many 
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’’equally likely" long-term traces of monthly streamflow or similar 
smoothly varying responses (Fiering, 1967) . 

Efficient management of water runoff from a watershed requires 
that hydrologic systems be described by dynamic models with suffi- 
cient accuracy. Since many controlling factors of hydrologic 
systems such as weather conditions, soil moisture variations, are 
known very little, or unknown, deterministic models do not at 
present offer satisfactory results. However, in this case, with 
certain reasonable assumptions, one can use random variables to 
approximate the stochastic nature of the system, and, then, can 
analyze it if one can track these "random" variables with time. 

For dynamic systems which are well characterized by finite- 
order ordinary differential equations (differential systems) with 
additive noise terms, when the analysis in the time-domain is to 
be preferred, the use of the so-called state— space approach will 
offer a great deal of convenience conceptually, notationally , and, 
sometimes, analytically. The use of state— space concepts and modern 
control and estimation techniques for the analysis of lumped-para- 
meter response models of hydrologic systems, when input and output 
data are corrupted by additive noises, will be discussed further in 
subsequent chapters. 

Identification Problems in Hydrology 

Zadeh (1962) defines the problem as 

"Identification is the determination, on the basis of input and 
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output of a system within a class of systems (models) , to which 

the system under test is equivalent (in some sense)." 

Thus, techniques for system identification must be based on 
data. For deterministic models, errors in data are reflected in 
the identification results and in errors in predicted outputs from 
an incorrect model of the system. An empirical study of the 
response of a simulation model of the rainfall-runoff process to 
input and output errors has been made by Dawdy, Lichty, and Bergmann 
(1972) . For their model, they found that the errors of streamflow 
estimates are approximately linearly related to errors of rainfall 
input data. 

Errors in rainfall-runoff data can result from the following 
sources (Rao and Delleur, 1971) : (a) errors in reading stage 

hydrographs, (b) errors inherent in the rating table, (c) errors 
associated with the method used in base-flow separation, (d) errors 
resulting from an inadequate network of precipitation stations and 
(e) errors in the method of determining the rainfall excess. While 
the errors associated with reading the stage hydrograph can be 
estimated, at present no accurate assessment can be made of the other 
errors involved. Since precipitation changes rapidly in time and 
space over a watershed and a network of precipitation stations may 
not be adequate, point rainfall measurements are unlikely to be 
representative of the actual rainfall on the watershed. On the 
contrary, the runoff is eventually collected at a single point, the 
mouth of the watershed, and the discharge data are usually of very 
good quality. It could, therefore, be reasonably assumed that the 
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runoff is reasonably noise-free compared to the rainfall. Thus, 
techniques that can use the observed outputs to compensate the 
disturbance in the inputs (i.e., control the inputs) and at the 
same time adaptively learn the characteristics of the unknown noise 
in order to get better system identification are needed in the 
study of the rainfall-runoff process. These techniques can easily 
be derived from modern control and estimation theory, based on 
state-space concepts. 

For large scale hydrologic systems it is not often possible to 
specify the a priori structure or functional form of the model 
(Kisiel, 1969). Structural information may be none, partial or 
complete. In the complete-information case, the system identification 
task reduces to a parameter estimation problem. The study reported in 
this thesis deals only with the applications of the state- variable 
approach in modem control and estimation theory to the identification 
of unknown parameters of nonlinear lumped— parameter response models of 
hydrologic systems subject to noisy input-output data. 

Advantages of the State- Variable Approach 

Modem control theory is based on the state-space concept. The 
idea of state as a basic concept in the representation of systems 
was first introduced in 1936 by A. M. Turing (Tou, 1964). Later, 
the concept was employed by C- E. Shannon in his basic work on 
information theory. The application of the state-space concept in 
the control field was initiated in the forties by the Russian scien- 
tists M, A. Aizerman, A. A. Fel'dbaum, A* M. Letov, A. I. Lur^e 
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and others. In the United States the introduction of the concept 
of state and related techniques into the optimum design of linear 
as well as nonlinear systems is due primarily to R. Bellman, The 
basic work of R.E. Kalman in estimation and control theory, and 
the extension of his work by others, played an important role in 
the advancement of modern control theory. 

The state of a dynamic system is defined (Ogata, 1970) as 
the smallest set of variables, called state-variables , such that 
the knowledge of these variables at t = t^ together with the 
input for t ^ t^ completely determines the behavior of the system 
for any time t ^ t^. The state-space representation of a system 
is not unique. In the design of optimum control systems it is 
extremely desirable that all the state variables be accessible 
for measurement and observation. For a linear system in the linear 
filtering problem, Athans (1967) showed that the choice of the 
state-variables is not crucial since one can obtain the estimates 
of another set of state variables using a simple linear nonsingular 
transformation. The same linear transformation also links the error 
covariance matrices in the two models. 

The advantages of the state-space concept over the conventional 
transfer function approach can be listed as follows (Ogata, 1970): 

(1) The state-variable formulation is natural and convenient 
for computer solutions . 

(2) The state- variable approach allows a unified representa- 
tion of digital systems with various types of sampling schemes. 
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(3) The state-variable method allows a unified representation 
of single variable and multi-variable systems. 

(4) The state-variable method can be applied to certain types 
of nonlinear and time-varying systems. 

Control Systems Terminology 

A system to be controlled, called a plant , has a set of 
outputs represented by the vector Y and a set of inputs represented 
by the vector U. A priori information about the plant may also be 
available, which usually is the desired output and represented by 
the vector R. 

Definition 1 ; An open-loop control system is one in which the 
control action is independent of the output. In this case, U is 
obtained as an operation on R, and the operator is known as an 
op^en-loop controller . The block diagram of an open-loop control 
system is illustrated in Figure 1-1. 


R 

open-loop 

-- u , 

PLANT 

Y 

REFERENCE INPUT * 

controller 

CONTROLLED 

INPUT 

CONTROLLED ^ 
OUTPUT 


Figure 1-1: Block diagram of an open-loop control system. 

Two outstanding features of open— loop control systems are: 

(1) Their ability to perform accurately is determined by their 
calibration. To calibrate means to establish or re-establish the 
input-output relation to obtain a desired system accuracy. 

(2) They are not generally troubled with problems of 


instability. 



9 


Definition 2 : A closed-loop control system is one in which the 

control action is somehow dependent on the output. Closed-loop 
control systems are more commonly called feedback control systems . 
In this case, U is an operation on R and Y, the operator is called 
a feedback controller . The block diagram of a feedback control 
system is illustrated in Figure 1-2. 



Figure 1-2: Block diagram of a feedback control system. 

When the summing point is a subtracter, i.e., e = R - b, one 
has negative feedback . When it is an adder, i.e., e = R + b, one 
has positive feedback . 

The most important features the presence of feedback imparts 
to a system are the following; 

(1) Increased accuracy. For example, the ability to faithfully 
reproduce the input . 

(2) Reduced sensitivity of the ratio of output to input to 
variations in system characteristics. 

(3) Reduced effects of nonlinearities and distortion. 

(4) Increased bandwidth. The banwidth of a system is that 
range of frequencies (of the input) over which the system will 
respond satisfactorily. 
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(5) Stabilized effect to an unstable system, i.e., without 
the addition of feedback, the unavoidable uncertainties in initial 
conditions and the inaccuracies in the model that would be used for 
determining an open- loop control would render such a system useless. 

If everything about the environment and process is known 
a priori, the design of the control law is straightforward and can 
be accomplished by means of proven techniques. On the other hand, 
if the environment or process is poorly defined, more advanced and 
sometimes less-proven techniques must be used to design the law. 

In the latter situation, control specialists have devised adaptive 
control systems and learning control systems. 

Definition 3 : An adaptive control system is one which is provided 

with: (1) a means of continuously monitoring its own performance 

relative to desired performance conditions, and (2) a means of 
modifying a part of its control law, by closed-loop action, so as 
to approach these conditions (Cooper and Gibson, 1960) . 

A comprehensive survey of adaptive -control systems was 
presented by Aseltine et al. (1958) . The definition of adaptivity 
and characteristics of adaptive control systems have also been 
treated by Zadeh (1963) , Braun (1959) , Donalson and Klshi (1965) , 
Eveleigh (1967) and Davies (1970) . 

Definition 4 ; A learning control system is an improved adaptive 
system which can memorize the optimal control function once 
established through adaptation and can immediately execute optimal 
control without adaptive search when a once experienced situation 
takes place again (Tamura, 1971). 
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A learning control system, therefore, should have memory 
facilities to store pairs of experienced situations and the results 
of adaptation. Moreover, it should have the capability to relate a 
certain control function with the present situation. A representa- 
tive learning control system can be illustrated as in Figure 1-3. 



input o'itput 


Figure 1-3: A representative learning control system. 

The state-of-the-art of learning control theory and applica- 
tions have been given by Gibson (1963) , Fu et al, (1963) , Mosteller 
(1963) , Sklansky (1966) , and Mendel and Zapalac (1968) . In their 
studies, various learning control models have been discussed and an 
extensive bibliography on the subject was presented- 

Objectives 

The main objective of the study reported in this report is 
the introduction of two approaches, namely adaptive control and 
sequential non-linear estimation, for the identification of the state 
and unknown parameters of a nonlinear hydrologic system response model. 
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A state-space implementation of these techniques was used 
for the analysis of the Prasad equation for the rainfall- 
runoff process. The algorithms developed are a Kalman filtering 
scheme with adaptive estimation of the error-covariance matrices 
and secondly, an interated extended Kalman filter. Two rather 
general computer programs were developed during the investiga- 
tion and are discussed in detail. 


Report Outline 

Following is the outline of the contents of each chapter 
in this report. Chapter 2 is devoted to a review of parameter 
identification techniques used in hydrology in the past. These 
techniques range from optimum search, recursive least-squares, 
to a few more sophisticated optimization methods using sensi- 
tivity analysis and quasi-linearization to identify unknown 
time-invariant parameters in a nonlinear model. Chapter 3 
presents the formulation of the identification problem, first, 
as a regulator problem in stochastic control theory, and then, 
as an optimum sequential estimation problem in a noisy 
situation; both methods use state-space concepts. The results 
of this chapter are two adaptive identification algorithms which 
can track unknown parameters in a nonlinear model with noisy 
observations. Finally, the implementation of the two proposed 
approaches to the study of the rainfall-runoff processes in a 
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selected watershed and a discussion of the results are presented 
in Chapter 4, Chapter 5 presents a summary of the study and 
some recommendations for further investigation. The derivation 
of the necessary equations mentioned in various chapters and 
a description of the computer programs that were written to 
implement the two approaches are presented in the attached 
appendices . 
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CHAPTER 2 

REVIEW OF PARAMETER IDENTIFICATION TECHNIQUES IN HYDROLOGY 

When the structure of the model of a hydrologic system is 
known, the system identification problem becomes the problem of 
estimation of unknown parameters in the model. The following devel- 
opment summarizes the various approaches proposed in the past to 
solve the parameter identification problem in hydrology. 

Optimum Search Techniques 

Search techniques are very useful for engineers and hydrologists 
to solve optimization problems when it is impossible or impractical 
to solve them directly by analytical optimization techniques. The 
only requirements are that a value of the function can be determined 
for any given set of variables and that, when a global extremum is 
sought, the function has no unbounded or multiple peaks (well-behaved 
functions) . When constraint equations on the variables are 
associated with a given problem, the objective or cost function may 
be augmented by penalty functions such that the extremum of the 
augmented but otherwise unconstrained functions converge to the 
cont rained extremum of the cost function in the limit, and the usual 
search techniques may be applied with little modifications. This 
very useful penalty function concept was first introduced by Courant 
(1943) and later modified by Carroll (1961) and by Goldstein and 
Kripke (1964) . The penalty argument has the defect that it may 
yield fictitious solutions when the problem is ill-posed. More 
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discussion on this subject has been given by Beltrami (1970) . A 
great number of search techniques have already been proposed in the 
past. One can easily find them in various textbooks or reports on 
optimization theory and control engineering (Davidson, 1959; Norris, 
1960; Wilde, 1964; Leon, 1966; Pierre, 1969), or in various technical 
journals dealing with computational methods in optimization (Brooks, 
1959; Rosenbrock, 1960; Powell, 1964; Shah et al., 1964; Fletcher, 
1965; Young, 1965). 

Search techniques can be grouped into two broad categories: 
deterministic search and random search. Techniques belonging to the 
latter category are superior in solving optimization problems of 
complex nonlinear hydrologic systems, such as rainfall-runoff 
processes, where discontinuities of the first derivatives and noise 
in the system can cause deterministic algorithms to become inefficient 
or to fail. Practical algorithms for random search have been proposed 
or discussed by Brooks (1958) , Hooke (1958) , Gurin and Rastrigin 
(1965), Gurin (1966), Schumer and Steiglitz (1968), Zakharev (1969) 
and Hill (1969) . Good survey papers on random search methods were 
also given by Kamopp (1963) and White (1972) . For general comparison 
purposes, search techniques can be divided into two classes: those 

techniques which utilize derivatives of the performance measure, and 
those techniques which do not. In general, the best sequential 
search techniques are more efficient than the best nonsequential 
ones, and the best sequential methods which utilize the gradient are 
more efficient than those which do not (Pierre, 1969). Search 
techniques, especially gradient methods, were used very often by 
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hydrologists for fitting parameters in a conceptual model of 
catchment hydrology. Ibbitt (1970) has tested eight deterministic 
optimum search methods and one random search method for application 
to hydrologic models. Among deterministic optimum search methods, he 
found that Rosenbrock's method was the best for the following 
reasons: (a) it was the most efficient among the methods tested; 

(b) it had an extremely flexible constraint technique; (c) its 
demands for computer storage were small; and (d) it could be applied 
to any type of objective function. 

Least-Squares Procedure 

Identification techniques based on least— squares procedure are 
applicable to both linear and nonlinear systems. The method of 
least-squares was initiated by Karl Friedrich Gauss in 1795 but 
detailed description of this method was not published by Gauss until 
1809, in his book Theoria Motus Corporum Coelestium . Some basic 
ideas of Gauss about the method of least-squares are: 

(1) minimum number of observations are required for the determi- 
nation of the unknown parameters; 

(2) model equations must be exact descriptions of real systems; 

(3) observation errors are unknown; and 

(4) the estimates of the unknown parameters must satisfy the 
observations in the most accurate manner possible. 

The best estimates of the unknown parameters are defined as the set 
of values that minimizes the sum of the squares of the observation 
residuals. Based on the least-squares concept, a recursive least- 
squares approximation algorithm can be formulated as follows (Duong, 1970) . 
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Given a time-invariant nonlinear hydrologic model 

= f(X,P) 

where y is the computed output (scalar) from the model; 


(2-1) 


X is the set of state variables of the system; and 

P is the set of unknown parameters to be identified; 

let y be the observed output from the system. The output residual 
o 

at sampled time t^ is defined to be 

~ (y ) . “ (y ) - (2-2) 

•^1 -^o X c X 

The best estimates of all unknown parameters of the model 

are computed at the same time by proceeding as follows; 

After guessing initial values, P°, for the unknown parameters, 

one writes Ay , in terms of the corrections Ap , , related to the 
1 J 

parameters p^ , as 


- 


I 

j=l 




( 2 - 3 ) 


where n is the total number of unknown parameters to be identified. 
In the least-squares regression method, the best corrections Ap^ 
to the a priori known values of the parameters will minimize the 
function 


J = (AY - AAP*) (AY - AAP*) 


( 2 - 4 ) 


where the superscript T denotes the transpose of a matrix, AY 
is the output residual vector whose elements are the Ay^, and A 
is the matrix of the partial derivatives of dimension N x n (N is 
the number of observations, N>n) , with elements 
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a. . 


J j j. 


1=1,2, ,..,N ; j=l,2,...,n (2-5) 


Minimization of the function (2-4) relative to AP* leads to the 

T 

unique solution, if A A is non-singular, 


* . T -1 T 

AP = (A A) A AY 


(2-6) 


and the best estimates of the unknown parameters p^ are given by 

5*f o * 

p = p'^ + AP (2-7) 

•ft 

The pj are then substituted back into equation (2-1) and the 

(y^)^ are computed again. New residuals are formed from equations 

(2-2) and new corrections to the p^ are computed from equation 

(2-6), which provides the next approximations to the parameters. 

One important assumption of least-squares regression technique is 

that the unknown parameters must be Independent. If this assumption 

is not satisfied, i.e., the "independent" variables are not truly 

independent, then the correction Ap^ will not be uniquely associated 

with p^ and convergence of the method cannot be insured. 

To account for the difference in accuracy that might exist 

between various measured outputs and possible relations between them, 
* 

the values p (best estimates of the p ) to minimize the function 
•J J 


J* = (AY - AAP*)^ W (AY - AAP*) 


( 2 - 8 ) 


are often sought, where W is an NxN weighting matrix and it is 
assumed to be symmetric and known. In general, W is chosen to be 
the inverse of the covariance matrix of the errors. If AW A is 
non— singular , the unique solution for the optimization problem 
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will then be given by 

AP* = (a”^W A)~^a'^WAY (2-9) 

The method of least-squares has been used extensively in 
hydrology during the last two decades for fitting parameters of 
parametric models of hydrologic systems. The recursive least- 
squares approach and the method of differential correction introduced 
to hydrologists by Snyder (1962) and later refined by Decoursey and 
Snyder (1969) to optimize hydrologic parameters are part of a more 
general theory on sensitivity analysis originated by Bode (1945) . 

The mathematical background of the theory of sensitivity analysis of 
dynamic systems was treated by Tomovic (1962) and can be found in 
many textbooks in control engineering (i.e., Perkins, 1972) and 
therefore is omitted here. Sensitivity analysis has been used by 
Vemuri et al. (1967, 1969) in the analysis of ground-water systems. 

Quasi-linearization 

In this section, we consider the possibility of solving a 
non-linear hydrologic system problem by first transforming it into 
a related linear system problem whose solution is then modified to 
obtain the desired solution. Only one such indirect computational 
method for solving system identification problems is treated here. 
This technique is known as quasi-linearization and has been used by 
hydrologists in the identification of a non-linear hydrologic 
system response (Labadie, 1968) and of unconfined aquifer parameters 
(Yeh and Tauxe, 1971). 
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The quasi-linearization technique, which is often referred 
to as a generalized Newton-Raphson technique, was originally 
presented by Bellman and Kalaba (1965) . Sage and Burt (1965) , Sage 
and Smith (1966) , Sage and Melsa (1971) and Graupe (1972) have 
examined the application of discrete and continuous quasi-lineari- 
zation to system identification problems. Following is the summary 
of the development of the continuous quasi-linearization technique. 
Consider a non-linear hydrologic system described by a vector 
differential equation of the form 

X = F(X,P,U) (2-10) 

where X is an n-dimensional state-vector; 

P is an r-dimensional unknown parameter- vector to be identified; 
and 

U is an m-dimensional input-vector. 

The n+r boundary conditions of the equation (2-10) are assumed to 
be linear and known, and the elements of P are stationary, i,e., 

P = 0 (2-11) 

The two equations (2—10) and (2— H) can be combined to have the 
form 

Z(t) = c[z(t),U(t)] (2-12) 

with boundary conditions : 


H(t)Z(t) 


B(t) 


(2-13) 
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where Z is the new (augmented) state vector defined by 

Z(t) = [x^(t) , . . . ,x^(t) ;pj^, . . . ,p^J (2-14) 

Expanding G^Z(t),U(t )] in a Taylor series about the i estimate 
of the state vector, Z^(t), the (i+1)^^ estimate of the state is 
then given by 

= G[z^(t),t] ^g[z (t),t] _ Z^(t)] 

izht) 


+ higher order terms. 


(2-15) 


Assuming that the initial estimate is close to the predicted value, 
and dropping the higher order terms in equation (2-15) , one obtains 

Z^'^^(t) = + G[z^(t),t] - 

^Z^(t) 


^c[z^(t) ,t. i 


^Z^(t) 


Z (t) 


This equation has the form: 


Z^'^^(t) = A^(t)Z^"^^(t) + v’^(t) 


i+1 


(2-16) 


where 


iz\t) 


(2-17) 


ht) =GfzV),t] - Z^t) 


(2-18) 
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The general solution of equation (2-16) has the form 




(2-19) 


where ^ (t,t^) is the fundamental solution of equation (2-16), 

given by 


' 1+1 

rl) 


f i- f- ) = /h 


i+1 


t- 1- 


( 2 - 20 ) 


with 




( 2 - 21 ) 


i+1 

and Q (t) is the particular solution of equation (2-16) , 
satisfying 

Q^'*'^(t) = A^(t)Q'‘^'^^(t) + V^(t) 

with ) = 0 

o 


( 2 - 22 ) 

(2-23) 


Substituting the general solution (2-19) into the boundary conditions 
(2-13) , one obtains 


H(t) [ 3>^‘^^(t,t^)Z^"''^(t^) + Q^'^^(t)] = B(t) 
H(t) $^"'’^(t,t^)Z^"‘'^(t^) = B(t) - H(t)Q^'*'^(t) 


which is of the form 


where 


A Z^'^^(t^) = B 

A = H(t)4.’-'''^(t,t ) 

o 


i+1 


B = B(t) - H(t)Q (t) 


(2-24) 

(2-25) 

(2-26) 

(2-27) 

(2-28) 
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i+1 

Thus, the solution for Z (t) has been reduced to a set of linear 
initial condition problems given by Eqns. (2-26) which are easily 
solved . 

The quasi-linearization technique for system identification 
often requires a good initial estimate of the states in order to 
converge. The computational effort in identification by this 
technique is considerable, and the approach is limited mainly to 
cases where only some states (not necessarily the same states) are 
accessible at different times (Graupe, 1972) . 

Summary 

In this Chapter, a review of various techniques for the identi- 
fication of unknown parameters of lumped-parameter models of hydrolo 
gic systems is presented. These techniques range from optimum 
search methods, least-squares procedure to more sophisticated optimi 
zation techniques using sensitivity analysis and quasi-linearization 
to identify unknown parameters in a non-linear model of a hydrologic 
process. Most of these techniques are suited for the analysis of 
linear or non-linear deterministic time- invariant systems; some of 
them can be used when the inputs are imbedded in noise; but none of 
them are good for the analysis of non-linear time-varying systems 
with noisy observations- For this case, adaptive learning control 
techniques and nonlinear filtering approaches, using state-space 
concepts, might be more suitable. Two such techniques will be 
presented in the next Chapter. 
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CHAPTER 3 

STATE-SPACE APPROACH FOR THE IDENTIFICATION OF NONLINEAR 
HYDROLOGIC SYSTEMS FROM NOISY OBSERVATIONS 

As already mentioned in Chapter 2, techniques used in the 
past for the determination of the instantaneous unit hydrograph and 
the identification of unknown parameters in a conceptual model of 
a hydrologic process are not adequate. The main reasons for this 
come from the fact that the input and output hydrologic data are 
imbedded in noise, the hydrologic processes are more or less non- 
linear, and the changing of the environmental conditions with time 
may affect the model output. In this chapter, a new approach using 
state-space concept is investigated, and techniques for optimal 
adaptive identification of unknown parameters and of the control 
inputs are presented. 

Problem Formulation 

A general nonlinear lumped-parameter model for a hydrologic 
system can be represented by the following vector differential 
equation 

♦ 

X(t) = f [x(t),U(t),t] + w(t) (3-1) 

where X(t) is the actual (nxl) state-vector of the system, U(t) is 
the noisy forcing input (i.e., rainfall), f(0 is the nonlinear 
mapping from E into E , and w(t) is an (nxl) noise term which is 
assumed to be a zero— mean white noise process with covariance 



25 


matrix 

E {w(t)w(x)^j = W 6(t - T) (3-2) 

where 6(») is the Dirac-delta function and superscript T indicates 
transpose . 

In most cases one cannot measure the state directly and 
precisely, but one can measure a vector Y which represents the 
output (i.e., runoff) of the system, and is related to the state 
by the following equation 

Y(t) = h[x(t)] +v(t) (3-3) 

ti X 

where h(.) is the nonlinear mapping from E into E , and v(t) is 
the observation noise which is again assumed to be a zero-mean white 
noise process with covariance matrix 

E Jv(t)v(x)^j = V 6(t - t) (3-4) 

The identification problem can now be stated as follows: 

Given the system model (3-1) , the observation equation (3-3) , 
the noise characteristics defined by Eqns. (3-2) and (3-4), and the 
set of noisy input-output pairs and ; from the initial 

conditions of the state and the estimation-error covariance matrix, 
find the best estimate of the state of the process under some 
optimality criterion (i.e., minimization of the mean-square-error of 
the estimate) . 

Two approaches for solving this problem are presented . In the 
first approach the given problem is treated as a regulator problem 
in stochastic control theory. In the second approach a non-linear 
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filter will be used to estimate the state and unknown parameters 
of the system under noisy observations. 

A Linear Stochastic Control Problem 


In order to apply linear stochastic control theory to the 

identification problem formulated above, one must linearize the 

given nonlinear system equations and the observation model to get 

the linearized model of the system. 

Let the nominal values of the state and the input be denoted 
* * 

by X (t) and U (t) respectively and the deviations of the actual 
state and the input from nominal values be represented by 


6x(t) = X(t) - X*(t) 
6u(t) = U(t) - U*(t) 


(3-5) 


The linearized equations of the system about the nominal values 
are then given by 

6x(t) = F(t) 6x(t) + G(t) 6u(t) + e(t) (3-6) 


where the elements of the matrices F(t) and G(t) are partial 
derivatives of f with respect to X(t) and to U(t) evaluated about 
the nominal values. 


F(t) 

G(t) 


^f(x.u.t) 

i X 


^f(x,u,t) 


I * * 

j x=x ,u=u 


I * * 

J x=x ,u=u 


(3-7) 


and e(t) denotes random disturbances which are used to model such 
effects as unknown dynamics and truncation errors. It is assumed 
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that e(t) is a zero-mean white noise process with covariance matrix 

E I e(t)e(x)^| = N6(t - t) (3-8) 

From now on, the vectors 6x(t) and 6u(t) will be referred to as 
the state and the control input, respectively, for the linear 
perturbation model. 

For calculations with digital computers, Eqn.(3-6) has to 

be converted into a difference equation. To that end one needs to 

have a state transition matrix <i'(t .,,t ) which satisfies the 

n+1 n 

homogeneous part of the differential equation (3-6); i.e.,4>(t ,,,t ) 

n+1 n 

is defined by 




(3-9) 


with the initial condition 


4>(t ,t ) = I for all t 

n n n 


(3-10) 


where I denotes the unit matrix having the same dimension as the 
matrix F. 

The difference equation of the system, derived from the differ- 
ential equation (3-6) for the time interval (t ,t .,), will have the 

n n+1 

following form 




t )6x(t ) + j 


^(t^+l> T)G(T)6u(t^)dT 



^(tn^l»T)e(T )dr 
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or, in abbreviated notations. 


6x - = ^ 1 6x + , 1 'Su + w 

n+1 n+l,n n n+l,n n n+1 


(3-11) 


The covariance matrix of the process noise w then becomes 

n 


E {w w } = 

n m 


ft 


n+1 


t 

J n 




nm 


= Q <5 
^n+1 nm 

where 6 is the Kronecker delta, 
nm 


(3-12) 


Similarly, linearizing the observation equation (3-3) about 
the nominal states and converting the result into a difference 
equation, one obtains: 


6y = H 6x + V 
n n n n 


(3-13) 


where 6y^ represents the deviation of the measured output from the 
nominal value, i.e.. 


6y = Y - Y (3-14) 

n n n ^ ' 

the elements of the mapping matrix H are partial derivatives of 

n 

h with respect to X evaluated at the nominal values: 


H 

n 


ih^(X) 



J x=x* 


(3-15) 


H is often called the observation matrix; and v is the observation 

n 

noise which is again assumed to be a zero— mean white noise process 


with covariance matrix 


T 

{v^v^ } 
n m 


R 6 
n nm 


(3-16) 
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In general, the Q and R are considered to be positive definite 


n 


and assumed to be given in advance. For complex systems, might 
be very hard to get. In this case, technique for adaptive estimation 
of is very useful and will be discussed later in subsequent 
sections . 

The control vector 6u in this optimal stochastic control 
problem is to be selected so that the performance index 

( N 


E = E 


N T T 

a , 6x^ + 6 u . , B. ^ 6u. 

jti j J ^ j 


(3-17) 


is minimized; where N is the total number of measurements made 
during the identification period, and the matrices A, and B. are 
arbitrary non-negative definite, symmetric weighting matrices. An 
approximate choice of these matrices must be made to obtain good 
results for the given problem. A choice that often turns out to be 

quite reasonable is (Bryson and Ho, 1969) 

_1 ^ 

A^ = n(t^-t^)xmaximum acceptable value of diag 6x(t)(Sx (t) ; 

—1 'p 

= m(t^-t^)xmaximum acceptable value of diag 6u(t)6u (t); 

where n and m are dimensions of the state-vector and the control 
input vector, respectively. 


The Optimum Controller 

In the above formulation, the equations (3-11) , (3-12) , (3-13) and 
(3-16) represent a linear time-varying system with Gaussian statistics, 
one can separate the solution of this optimum stochastic control 
problem into a deterministic optimum control problem and an optimum 
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estimation problem (Sorenson, 1968). The proof of the Separation 
Principle is given in Appendix A. The optimum stochastic control 
law is described by 




N-k \-k,N-k-l '^\-k-l 


(3-18) 


where A,, , is obtained from the solution of the deterministic 
N -k 

control problem, and is the optimum estimate of the state 

obtained by solving the optimum estimation problem. The 
are given by the following recursive formulas: 


^^-k N-k,N-k-l^N-k ^ N-k,N-k-l ^N-k-1 ] 


N-k,N-k-l^-k 

(3-19) 


^N-k ^N-k+l,N-k^N-k+l^N-k+l,N-k \-k 


(3-20) 


f I 

^N-k " ^N-k “ ^N-k ^N-k,N-k-l \-k (3-21) 

k = 0,1,2, . . . ,N-1. 

To start the iteration process, one can assume that ~ 0 

then proceed the calculations backwards in time to the initial 

control time t . The L„ , and , are often called the control 
o N-k IJ-k 

cost matrix and the control gain matrix , respectively. 

Identification of State— Variables 

The optimum feedback control law defined by Eqn. (3-18) 
depends on the optimum estimate of the state at each stage of the 
process under study. This optimum estimate can be obtained through 
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the well-known Kalman filter (Kalman, 1960) which is the best 

linear minimum variance estimator since the estimate 6x , defined 

n 

as the conditional expectation of given the available data 

set Z at time t , i.e., 6x =E {6x Iz }, is chosen to minimize 
n n n n * n 

the mean-square-error 

E{(6x -<5x)^(6x -6x)} = trace E{(6x - 6x )(6x - 6x 

nnnn n^^nn 

(3-22) 

The derivation of the standard Kalman filter is omitted here 
since one can easily find it in any textbook on estimation theory 
(Sorenson, 1966; Sage and Melsa, 1971) . Following is a summary of 
the filter algorithms for discrete cases where the state equation 
and observation model are given in the form of Eqns . (3-11) and 
(3-13), with noise covariance matrices defined by Eqns. (3-12) and 


(3-16) , 



one-stage 

prediction 

6x' / = 1 ^ J .1 

n+l/n n+l,n n n+l,n n 

(3-23) 

a priori 
variance 

^n+l/n ^n+1 ,n^n^n+l,n ^ ^n+1 

(3-24) 

Kalman 

gain 

^n+1 ^ ^n+l/n\+l[\+l^n+l/n^n+l ‘‘‘ ^n+l] 

(3-25) 

Filter 

algorithm 

'^''n+l " ^^n+l/n ^n+l['^^n+l ‘ “n+l'^'"n+l/n ] 

(3-26) 

a posteriori 

p — — K H P ' 

variance n+1 n+1 n+1^ n+l/n 

(3-27) 
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The initial conditions fix and P must be given, also the exact 

o o o » 

values of 0^ and must be known. After each observation, the 
optimum estimate of the state X is then given by 

X = X* + fix (3-28) 

The most important properties of a Kalman filter can be 
summarized as follows: 

(1) The filter estimates are all variables of the state vector 
in the least-mean-square-error sense. 

(2) The estimation is based upon statistical data of all error 
sources and is completely carried out in the time domain. 

(3) The filter formulae satisfy minimum variance criteria for 
all problem parameters. 

(4) The formulae implemented are recursive. This means that 
the optimum estimate at the present time can be computed from the 
previous estimate and the current observation without recourse to 
earlier estimates or observations. 

(5) The recursive formulae are well suited to digital computers. 

The use of the Kalman filter to estimate the unknown parameters 
in a lumped-parameter model of a hydrologic system is discussed in 
the next section, and the application of this approach to the analysis 
of rainfall-runoff is presented in Chapter 4. 

Parameter Identification 


Consider an unknown parameter a which is slowly varying in time. 
One could model its behavior satisfactorily by the random walk model 


a , = a + (w-) 
n+1 n an 


(3-29) 
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where (W~) is-a zero mean white noise sequence with covariance 
an 

matrix (Q-) . Equation (3-29) combined with Eqn. (3-11) gives an 
^ n 

augmented state equation of the form 


6x 



0 


6x 

+ 

V 

6u + 

w 

a 

n+1 

0 

1 


a 


,0. 

n 

w- 


(3-30) 


The observation model is also modified as 





V 

n 


(3-31) 


One could then apply the Kalman filtering algorithms mentioned above 

to estimate the unknown parameter along with values of the previous 

state variables of the system, using the augmented state equation 

and observation model defined by Eqns . (3-30) and (3-31). 

If the secular variation in a parameter a^ is rapid it can be 

decomposed into the product of a known, non-singular and rapidly 

time-varying transformation matrix, T , and an unknown but fixed or 

n 

only slowly varying parameter a . Thus, 

n 


and from Eqn. (3-29) : 


a = T a 
n n n 


a 

n 


$ a T + 
a n-1 


r (w-) . 

a a n-1 


(3-32) 


-1 p 

where <& = T T and 1 = T . One then uses Eqn. (3-32) to form an 

a n n— 1 an 

augmented state equation and new observation model, and proceeds 


as mentioned earlier. 
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Estimation of Error-Covariance Matrices 

In most practical situations , complete knowledge of the 
hydrologic process and measurement noise statistics is hard to get. 
The use of wrong a priori statistics in the design of a Kalman filter 
can lead to large estimation errors or even to a divergence of 
errors. The covariance matrix becomes unrealistically small and 
optimistic; the filter gain thus becomes small, and subsequent 
measurements are ignored. The state and its estimate then diverge, 
due to model errors in the filter. Analyses of error divergence 
in the Kalman filter have been performed by Heffes (1966) , Schlee et 
al. (1967), Nishimura (1967), Price (1968) and Fitzgerald (1971). 

Several approaches have been proposed for preventing filter 

. Horlick and Sward (1965) Investigated the technique of 
filter reset to keep the diagonal elements of the state covariance 
matrix above a specified value, Peschon and Larson (1965) proposed 
to include a random variable at the model input to account for any 
unrealistic assumptions made about the system model. Schmidt (1967) 
suggested two methods. In one method an estimate of the state is 
computed as a linear combination of the estimate given all prior 
data with the estimate given no prior data. Past information is 
degraded. The other method assumes a priori lower bounds on certain 
projections of the covariance matrix. Schmidt et al. (1968) also 
designed the modified Kalman filter equations by including an addi- 
tive gain matrix to the conventional gain matrix of the Kalman 
filter to prevent divergence of the filter resulting from unmodeled 
and computational errors. Recently, Tarn and Zaborszky (1970) 
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used a scalar weighting factor s for the observation error-variance 
and came up with a non-diverging filter which reduces to the 
standard Kalman filter for s = 1. 

The most interesting approach to prevent filter divergence 
is probably to cover model errors with noise and adaptively 
estimate the noise levels. The purpose of an adaptive filter is, 
then, to reduce or bound the noise by adapting the Kalman filter 
to the real data. A number of approaches to adaptive filtering has 
been proposed in the past few years. Good survey papers in this 
area were presented by Sage and Husa (1969), Weiss (1970) and 
Mehra (1972) . Following are some adaptive schemes for estimating 
the observation error covariance matrix R and the process error 
covariance matrix Q. 

1/Estimation of R 


One approach to reduce the effect of earlier measurement as 
new measurements are included is to replace the finite time 
averaging operation by exponentially weighted past averages . This 
technique, applied to the estimation of R^, is expressed as 


R 

n 


(1-W )R , + W R’ 
n n-1 n n 


(3-33) 


th 


where R^ is the predicted value of R at the n measurement: i.e.. 

n n * ’ 


R' = (6y - H 6x’ J(6y - H 5x ' , J 

n ^n n n/n-1^ ' -^n n n/n-1^ 


(3-34) 


This method has been applied successfully to satellite tracking 


problems, R is assumed to be constant during the sample period 
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For these problems, Tapley and Born (1971) proposed 

the following expression for the weighting coefficient W : 

n 

= (n-l)(n-2) ... (n-k)/n'^^^ (3-35) 

in which k is an integer. is zero for all nj<k, and as n-x», 
approaches 1/n. Hence, for n<k the value of R is not changed 
from the a priori value. The choice of k will depend upon how well 
is known, i.e., the more accurate R^ the larger k may be. 

In the case of a rainfall-runoff process, there are typically 
few data points obtained for one storm. By using the weighting co- 
efficient defined by Eqn. (3-35) one might not get the correct value 
for R. At the end of the observation period, the estimate R might 
still depend very much upon its initial guess. For this case, it 
is suitable to choose the weighting coefficient W as 


W 

n 


g 

l-d-a)"" 


(3-36) 


where cL ^ positive constant, normally less than 0.1 (Young, 1965). 
This choice of the weighting coefficient also assures that new data 
continues to have some effect on the estimate of R as long as the 
observation of the process is still going on. 


2/ Estimation of Q 


denote the predicted measurement residual, i.e.. 


A 

^n+l/n 


6v - H 5^' 

^n+1 n+1 n+l/n 


• I 
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one then has 


^ ^ ^n+l/n^n+l/n^ ^n+1 ^n+l/n ^n+1 ^n+1 


= Vl%+l.a"n*n+l,a<+l + VlVa+1 + \+l 


(3-38) 


Assume Q has the form 
n 


Q = q S S 
n ^ n n 

where q is a scalar; then, for the case of single output, Eqn. (3-38) 
can be written as 


E{ = H . , - P T , + 

n+l/n n+1 n+l,n n n+l,n n+1 


T T 

q H ,^S S H + R 
n+1 n n n+1 ^+1 


(3-39) 

To determine q, Jazwinski (1971) suggested that one finds the q 
value which satisfies the requirement 


max p (3_40) 

q>0 

or, when the probability density p is normal and has zero-mean, 


2 

^ n+l/n 


a+l/aJ 


(3-41) 


Equation (3-41) represents the consistency requirement for the 

estimates of q. From Eqns. (3-39) and (3-41), the following adaptive 

scheme is derived : 

f 2 T T 

^ n+l/n n+1 n+l,n n n+l,n n+1 n+1 , if positive 


^n+1 


T T 

H S H 
n+1 n n n+1 


(3-41) 


0 


otherwise 
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Since the estimate of q is based on only one residual, the 
estimator (3-41) will respond to (large) measurement noise samples 
as well as large residuals caused by model errors. One approach to 
remedy this defect is to apply the exponential smoothing technique 
described in the previous section, i.e., 


/I _ T.T \ A 

' "n+1" 


+ W 

n n+1 


An Adaptive Control Algorithm 


From the material in the above sections, the following 
algorithm can be set up to estimate the values of the state variables 
and unknown parameters of a time-varying system in a noisy situation. 

(1) Set initial conditions; 


"5^1/0 (usually set equal to zero) 



(2) Compute (off-line) the optimum control law by Eqns . (3-19), 
(3-20) and (3-21) from nominal values of the states. 

(3) Set n = 1. Start the sequential estimation process. 

(4) Adjust the observation error covariance matrix R bv 

n 

Eqns. (3-33) and (3-36), 

(5) Compute the filter gain by Eqn. (3-25). 

(6) Process the observation 6y^ by Eqn. (3-26) to get the best 

estimate 6x of the state, 
n 

(7) Replace by and return to step (4) to re-estimate 

R and 6:^ . 
n n 
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(8) Continue step (7) until one gets stable values for and 

(Sx . Hence, 
n 

X = X + (5x 
n n n 

(9) Compute the a posteriori error variance of the estimates 
by Eqn. (3-27) . 

(10) Adjust the model-error covariance matrix by Eqns. (3-41) 
and (3-42). 

(11) Compute the predicted state Eqn. (3-23) . 

I 

(12) Compute the predicted error variance of estimates 

by Eqn. (3-24) . 

(13) Set n = n+1 and return to step (4) . 

Three remarks can now be made about the above algorithm. 

First, during computation one must somehow design a technique 

to maintain the positive definiteness and symmetry of the matrix 

to avoid filter divergence. One way to maintain the desired property 

of the matrix P is to replace Eqn, (3-27) by an equivalent form 
n+J_ 

^n+l = - Vl + Vl Vl 

This equation is the sum of two symmetric positive definite matrices, 
and, when these are added, the sum will also be positive definite; 
therefore it is better conditioned for numerical computations than 
the previously mentioned expression (3-27) . 

Secondly, one could improve the above algorithm by correcting 
the nominal values after each new estimate of the state was obtained, 
i.e., one takes the new nominal value X as 
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X* = X* + 6x = ic (3-44) 

n n n n 

the new deviation of the state at stage n becomes 

= 0 (3-45) 

This process is often called trajectory rectification . As a conse- 
quence of relinearization, large initial estimation errors are not 
allowed to propagate through time, and, therefore, the linearity 
assumptions are less likely to be violated. 

Lastly, in order for the above algorithm to be called adaptive 
control, at each iteration step the new estimate of the state must 
somehow be used to alter the control law. However, the recursive 
equations used to compute the control history run backward in time. 

To update such a control law, the gains over the entire control 
interval would have to be recomputed each time a new estimate of the 
state was made available. The computational requirements of such a 
technique are enormous . A sub-optimum control law can be obtained 
as follows: first, the control law is computed off-line from nominal 
value of the state and the values of the matrix L are stored. Then, 
as each new estimate of the state is obtained, the control is updated 
by using new values of the matrices $ and r. The adaptive control 
algorithm can now be summarized as illustrated in Figure 3-1. 

When applying the above algorithm to the determination of 

unknown parameters in rainfall— runoff processes , the disturbance 

from rainfall data will be compensated by the optimum value of the 

control term, f ^ 6u , computed at each observation from the 
n'TX , n n 
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Figure 3-1: Block diagram of Adaptive Control Algorithm. 
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estimate of the state of the system and from Eqns. (3-18), (3-19), 
(3-20) and (3-21) . The adjustment for noise in rainfall data by this 
adaptive control approach is effective if the system is not highly non- 
linear, since the estimation of the best value for the state and the 
derivation of the optimum control law are all based on the validity 
of the linearized model. For this reason, another method for the iden- 
tification of unknown parameters in nonlinear lumped-parameter models 
of hydrologic systems, using nonlinear estimation techniques, is 
presented in the following section. The comparison of results of the 
two approaches, when applying to the identification of rainfall- 
runoff processes, is discussed in Chapter 4. 

A Nonlinear Filtering Problem 

A lumped-parameter model for the rainfall-runoff system is 
represented by Eqns, (3-1) , (3-2) , (3-3) and (3-4). Since rainfall data 
are noisy, one can approximate the correct rainfall input U(t) in 
the model (3-1) by 

U(t) = U*(t) + Wy*(t) (3-46) 

where II (t) represents the actual rainfall data, w^vsr(t) is the 
rainfall measurement noise which is assumed to be a zero-mean white 
noise process with covariance matrix 

T 

E {w^ft(t) Wy*(x) } = W^*6(t - t) 

The process equation is then written as 

X(t) = f [x(t),U*(t),t] + w(t) (3-47) 


■> 
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where w(t) is the combined noise term in the model which includes 
errors due to unknown dynamics and noise in input data. The form of 
w(t) depends on the form of w(t) and the relation of U (t) in the 
model equation. Assume that w(t) is also a zero-mean white noise 
process , having new covariance matrix defined by 

E {w(t) w(t)^} = W (t - t) (3-48) 

Since U (t) is a given scalar at time t, Eqn. (3-47) can be written 
simply as 

X(t) = f [X(t),t] + w(t) (3-49) 

The identification problem of the rainfall-runoff system then becomes 
the estimation of the state of a nonlinear system defined by the Eqns. 
(3-49) , (3-48) , (3-3) and (3-4), given the initial conditions X(0) and 
E{X(0)X(0)^} = P(0). 

Optimal estimation in the nonlinear case involves the solution 
of an infinite-dimensional process, as shown by Kushner (1967). Since 
the computational aspects of the truly optimum nonlinear filter are 
prohibitive, several approaches to sub-optimal filtering have been 
proposed in the past few years (Friedland and Bernstein, 1966; 

Schwartz and Stear, 1968; Athens et al.; 1968; Sage and Melsa, 1971). 
These algorithms can be roughly subdivided into the so-called first- 
order filters and higher-order filters with increasing complexity 
and computational requirements. Because in hydrology the estimate of 
the state of a system is usually not required to be highly accurate, 
only the extended (first-order) Kalman filter is considered here. 
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This filter is simple but effective and has been used very often in 
similar problems in the aerospace field. 

The Extended Kalman Filter 


This filter is the result of the relinearization procedure 

mentioned in the previous section. If, initially one linearizes the 

model equation (3-49) about X(t ) , then 

o 

6x(t ) = 0 

o 

The predicted deviation, given by 

6x'(t /t ) = ^>(t-,t ) 6x(t ) 

-L O to O 

is therefore equal to zero. 

Since one subsequently linearizes about X(t^) , 

6x(t^) = 0 

it follows that ~ ^ 

Thus, in general. 
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in view of relinearization, and using Eqn. (3-50), the correction 
to the estimate at an observation (Eqn. (3-26)) leads to 


\+l ^ n+1 ''' \+l f^n+1 ■ 


n+l/n 


>W)] »-52) 


Thus, the extended Kalman filter can be summarized in the following 
operations : 


A priori 
estimate 


x' 

n+l/n n I ^ 

J ^n 


f [x(t,t ),t] dt (3-53) 


P , = 4 . , (X ) P 4 , , (X ) + Q . , 

n+l/n n+l,n n n n+l,n n n+1 (3-54) 


A posteriori 
estimate 


= Vl/n-^ Vl<Vl/n>[Vr Vl^Cl/n^] 

(3-55) 

= - \+l('C+l/n> ' Cl/n 

•[^ - \+lK+l/n>f 

* Vl^Vl/n) Vl (3-56) 


Kalman 

gain 


^n+l^\+l/n^ 


^n+l/n \+l^\+l/n^ 


[\+l^\+l/n^ ^n+l/n \+l^\+l/n^ 




(3-57) 


The matrices and H are those of the linearized system defined by 


^ j-i (X ) 6x + w 
n+1 n+l,n n n n 


(3-58) 
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H (ic ) 6x 
n n 


n 


+ 


V 

n 


(3-59) 


The Iterated Extended Kalman Filter 


To improve the performance of the extended Kalman filter, one 
can use the technique derived by Denham and Pines (1966) to reduce 
the effect of measurement function (h) nonlinearity which occurs 
very often in hydrology when output data from a hydrologic system 
are imbedded in noise. This technique is a local iteration algorithm 
based on the relinearization about the new estimate. 


Consider the estimator, Eqn, (3-55), in the extended Kalman 
filter. It was obtained by evaluating the correction to the estimate 

'^n+l = K+X/n] <3-60) 

* T 

about Then, processing the observation via Eqn, 

(3-55) one gets Vi- Assuming that is closer to the true state 

t 

than one then would expect to get a better result by re- 

iincarizing the system equation about and recomputing the 

estimate. Thus, the iterated extended Kalman filter consists of Eqns. 
(3“53) - (3-57) with Eqn. (3-55) replaced by 


y(i+l) _ 
n+1 


= X 


n+l/n 


+ K -(x^^-h[Y - h 

n+1 n+1 


n+l^n+1' L n+1 


H 

n+1 n+1 n+l/n n+1 


>] 


i = 1,2, ,i (3-61) 


, ^(1) I 

where “ ^n+l/n’ local iteration terminates when there is 

no significant difference between consecutive iterations. The covari- 
ance matrix in Eqn, (3-56) is then computed based on the last estimate 
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of X 


Combining the sequential estimation of error-covariance matrices 
mentioned in the previous section with the above iteration for the 
improvement of the estimate of the state, the following algorithm is 
formed as illustrated in Figure 3-2. 

The local iteration process mentioned in this algorithm is 

designed for measurement nonlinearities and does not improve the 

previous nominal value chosen for the state on the Interval I * 

To include the nominal value in the iteration loop, one needs to 

smooth back the new estimate at t , , to t to get an improved nominal 

n+i n 

value for prediction to The linearized smoother is given by 

Jazwinsky (1970) as 


n/n+1 n n x L n+1 n+l/nJ 


with [Cl/nl 


( 3 - 62 ) 


or, using the smoothed estimate as » 


^i+1 = \ - Vl/n] 


(3-63) 


, A A f 1 ) t 

The iteration starts with C, = X and : = X , - , . 

1 n n+1 n+l/n 


This iterated linear filter-smoother, as named by Jazwinsky 
(1970) , was apparently first derived by Wishner et al. (1968) in a 
different way, and was called by these authors a "single state 
iteration filter". Although the performance of this algorithm was 
shown (Wishner et al., 1968) to be better than the iterated extended 



48 



Figure 3-2: Block diagram of the Iterated Extended Kalman Filter. 
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Kalman filter, the amount of computer time required is also two or 
three times greater. Hence this technique is suitable for those cases 
when one has only a small set of observation data and wishes to come 
up with acceptable estimates for the state of the system. For large 
data sets, the other algorithms can, hopefully, give the desired 
values for the estimates after processing all the data with a reason- 
ably small computer time. 

Summary 

In this Chapter, two approaches are presented for the estimation 
of the state and unknown parameters in a nonlinear lumped-parameter 
model of a hydrologic system. The optimum linear stochastic control 
approach is suitable for those cases when the error-covariance matrix 
of the input disturbances is unknown. The nonlinear estimation approach 
is simpler and more powerful, but requires knowledge about input noise 
characteristics. Subsequent chapters are devoted to the implementation 
of these techniques to the study of a particular model of the rainfall 
-runoff system. 
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CHAPTER 4 

IMPLEMENTATION AND RESULTS 


Nonlinear Lumped-Parameter Models for Rainfall-Runoff 

The nonlinearity of the rainfall-runoff relationship has been 
of concern only in the last decade; however, the concept of non- 
linearity and its methods of analysis are still very limited (Chow, 
1967) , Following are two nonlinear lumped-parameter models that have 
been used quite often by hydrologists in the past, namely the 
Kulandaiswamy model and the Prasad model. 


The Kulandaiswamy Model 


Direct runoff may be considered as the result of the transform- 
ation of rainfall excess by a basin system. The physical process of 
this transformation is very complex, depending mainly upon the storage 
effects in the basin. Kulandaiswamy (1964) derived the following 
general expression for the storage 


S = 


- a„(Q,R) ^ + f b^(Q.R) ^ (4-1) 


n=( 


n 


dt“ tsm 


dt 


m 


where S is the storage, t is the time, N and M are integers, and 
3-^(Q»R) end b^(Q,R) are parametric functions of the direct runoff Q 
the excess rainfall R. To apply Eqn. (4—1) to the study of the 
rainfall-runoff process in a particular watershed, the values of N 
and M must be determined. Both Q(t) and R(t) are available in the 
form of curves and differentiation has to be done by numerical 
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approximation techniques. Taking into consideration the nature of 
the curves representing Q(t) and R(t) and the magnitude of error 
likely to be introduced by numerical differentiation, the values of 
N = 1 and M = 0 have been adopted in Kulandai swamy ' s study. Eqn. 

(4-1) reduces to 

S = a^(Q,R) Q + a^(Q,R) ^ + b^(Q,R) R , (4-2) 

Plots of a , a^ and b versus Q , the peak discharge, for Willscreek 
o 1 o p 

Basin are illustrated in Figure 4-1. Kulandaiswamy found that a^ and 
b vary from storm to storm, but do not show any well defined trend 
in the variations; hence, he took these two parameters as constants 
(Kulandaiswamy and Subramanian, 1967) . The storage equation can now 
be written as 

S = a (Q) Q + a, ^ + b R . (4-3) 

o 1 dt o 

With the continuity equation 

II = R(t) - Q(t) , (4-4) 

the rainfall-runoff process can be represented by the following 
differential equation 

+ A(Q) + Q = (4-5) 

dt 

da 

where A(Q) = a + Q — ~ 

o dQ 

A plot of Q versus A(Q) was made for various basins and two types of 
regions could be differentiated. The system equations for these 




0.15 0.20 0.25 


QP, Inches/Hour 


Figure 4 1: Plots of a^, and vs. QP for Willscreek Basin 




5'J 


regions are 


(1) Non-linear region: 


a, ^ + (c^ + mQ) + Q 


dt^ 


o dt 


(4-6) 


(2) Linear region: 


2 


dt 


2 


dt 


+ Q = 


o dt 


(4-7) 


The general nonlinear storage equation (4-1) proposed by 

Kulandaiswamy has been accepted by many hydrologists in the simulation 

of the rainfall-runoff process by lumped-parameter response models, 

but the approach used in the determination of the model parameters 

has also been criticized (Eagleson, 1967) - Kulandaiswamy used 

characteristics of the surface runoff hydrograph at peak discharge 

(dQ/dt = 0) , on the falling limb (R = 0) , and on the rising limb up 

to the end of rainfall excess to get various plots of a , a- and b 

o 1 o 

vs. Q and Q vs. A(Q) ; then from these plots the values of a., c- , m, 

P J. i 

b and c^ were determined. The evaluation of a from a single 
o 2 o 

discharge ( the peak discharge ) and a^, b^ from a portion of the 
surface runoff hydrograph should be replaced by some other means 
that can evaluate the model coefficients over the full range of 
observed discharges. 


The Prasad Model 


A simplification of the above model by retaining only two terms 
of the general nonlinear storage equation was proposed by Prasad 
(1967); in this case the storage equation is 
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S = + K2§ (4-8) 

in which may be a complicated function of several variables 
affecting the wedge-storage as well as the storage-discharge relation- 
ship. In his study, Prasad assumed that K^, and N are constant 
for a particular hydrograph. Using the continuity equation (4-4), one 
gets the following differential equation for the rainfall-runoff 
process 

Kz + KiN ^ + Q - R . (4-9) 

Comparing the Prasad model for nonlinear storage (Eqn. (4-8)) with 

the Kulandaiswamy model defined by Eqn. (4-3), one can recognize 

that a^(Q) and a^^ have been taken as K^Q ^ and respectively, 

and b =0. 
o 

In the Prasad model, the time-invariant coefficients 
and N were evaluated by a trial-and-error method which is computation- 
ally inefficient and requires the knowledge of the initial conditions 
with sufficient accuracy. These coefficients were later computed by 
Labadie (1968) using quasi-linearization technique which has two 
main inherent weaknesses: (1) Initial approximations must be within, 
or at least close to, the convex region surrounding the optimal 
solution, or convergence is not attained. (2) If convergence does 
not result for a particular set of initial approximations, it is not 
possible to determine systematically a better set of initial approxi- 
mations from these results. 

^11 the above approaches for estimating the model coefficients 
are suitable only for deterministic models and not suitable for the 
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analysis of real input-output data where the values to be used in 
the model are imbedded in partially known or unknown noise. The two 
methods proposed in Chapter 3 are very useful in solving parameter 
identification problems in this case. 

Since the Prasad model for rainfall-runoff is typically non- 
linear and the data set made available to the author was related to 
Prasad's work, only the Prasad model will be used in the investigation 
of the proposed identification schemes' performances. 


Reformulation of the Prasad Model in State-Space 
Eqn. (4-9) can be written as 


dt^ 


-(| )KiN § - (i )Q + (| )R 


Using the transformation 


X. 


Q 

Q 

^1 

1/K, 


= N 


(4-10) 


(4-11) 


and the assumption that the model coefficients are t ime- invar i an t , 
Eqn. (4-10) can be written in the following form 


* 


• 




,1 


Xc-1 2 

^2 


-WsV X2 + X^(R-xp 

^3 


0 



0 

X- 


0 

5j 




(4-12) 
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or, in abbreviated notation, 

X(t) = f (x(t),R(t)] . (4-13) 

Eqn. (4-13) is the model equation in state-space. Let Y(t) denote the 
measured runoff which is embedded in noise, one then has 

+ v(t) , (4-14) 

or, in abbreviated notation, 

Y(t) = h [x(t)] + v(t) (4-15) 

where v(t) represents the noise tern. 

Eqns. (4-12) , (4-13) , (4-14) and (4-15) are the basis for 
further development. 

Computation of the System State— Transition Matrix 

The crucial problem in applying the proposed estimation schemes 
to continuous systems with discrete measurements is the evaluation of 
the state-transition matrix. 

It is known that can be obtained from the coefficient 

o 

matrix F(t) by the differential equation 

*t,t = 

o o 


Y(t) =[ 10000 ] 


with the initial condition 
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$ ^ = I 

t ,t 
o o 


(4-17) 


In the general case an explicit closed form solution of Eqn. 
[4-16) is not possible, but an infinite series can be derived which 
:onverges uniformly in t for every matrix F(t). This solution is 
;iven as follows, by the use of the Peano— Baker method of integration 
[Pipes, 1963) 




^ O y 


F(T2><iT2^'^l + 


A 


F(x^) 


J 


A A 

^ FCt^) ^ F(T^)dT^dT2dT^ + ... (4-18) 

t t 

o , o 


[f the matrix F(t) satisfies the commutativity condition, i.e,, 


F(t^) FCt^) = ^(^ 2 ) F(tj^) ^2 


:hen the state-transition matrix is given by (DeRusso, Roy and Close, 
L967) 


^ - exp 

t, t 

o 


F(T)dT 

t 

; ° 


I I + ( F(T)dT + i [ 

1 > V 

+ Y, t { F(T)dT + . 
* t 

O 


F(x)dT ] 


(4-19) 


Sxpansion of each term and rearrangement yield: 




2 3 

At • 2 At * * 3 • • 3 

. = J I + At.F +^(F + F)+^[F + -(FF + FF)+F] 

o 2! o o 3! o 2 o o 00 o 


} 


(4-20) 


jfhere F = F(t ) and At = t - t . 
0 0 o 
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The computation of Eqn. (4-20) contains only a sufficient 
number of terms so that additional terms are negligible by comparison 
with the partial sum to that point. However, there are some diffi- 
culties associated with attempting to use a truncated form of this 
expansion. These stem from the impracticality of obtaining the first- 
and higher-order derivatives of F(t), if second- or higher-order 
accuracy is required. In these cases, it is better to subdivide the 
interval (t - t^) and to consider F as a constant matrix during 
these partial Intervals. The partial transition matrix can now be 
computed by 


^ = I + F At,, + ^ + 

h’h-1 ^ 2 ! 



(4-21) 


where At^^ - (At)/N and N is the number of sub-intervals. The state- 

transition matrix $ over the whole interval (t - t ) is simply 

t,to o 

the product of the partial matrices ^ 

^i’^i-1 



$ 







(4-22) 


t = 


t > t 1 > t « > 
p p-1 p-2 


> t- > t 
1 o 


For fast time-varying systems, it is better to use small sub- 
intervals and consider only the linear term in the expansion (4-21) 
than use bigger sub-intervals and include higher-order terms in the 

computation of $ (Unger and Ott, 1970). 

o 

The formulation of the estimation schemes for estimating 
system parameters with noisy input-output data have been implemented 
for the study of the rainfall— runoff process. The data to be used 
were from the storm of April 10, 1953, on South Fork Vermilion River 
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Basin above Gatlin, Illinois. The following results are obtained 
from computer programs written in EXTENDED FORTRAN for use on the 
GDC 6400 digital computer at Colorado State University. 


Adaptive Control Approach 


The matrices F(t) and G(t) in the linearized expression of 
Eqn. (4-12) have the following forms 


F(t) 


where 


0 

^1 

0 

0 

0 

^1 = 

E„ = 


0 


X 1 

-XaX^X^Xi -X2X^X3X^ 


X,_1 


X5-I 


X.-l 


E. = - 


X3-I 


^2^3 Vl » 

G(t) = ^0 0 0 oj^ 


(4-23) 


(4-24) 


The error term e(t), used to model the truncation error, is 
chosen to be a zero— mean white noise process with covariance matrix 

0 


Q = 


.01 


0 


0 


0 


The observation matrix derived from Eqn. (4-14) has the following 


form 
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H = 



(4-25) 


the observation noise v(t) is also assumed to be a zero-mean white noise 
process with variance 

R = .0001 


R is taken to be much smaller than Q because the observed outputs 
of the rainfall-runoff system are relatively noise-free compared to the 
inputs; also, there are errors in the model equations due to the incom- 
plete knowledge of the nature of the system. 

The following initial conditions are assumed for starting the 
adaptive control algorithm: 

X*(0) = X*(0) = 0. 

X*(0) =10. , X^(0) = .1 , X*(0) = 1. 

=0, , i = 1, 2, 3, 4, 5 


P = 
o 





The initial estimation-error variances are chosen somewhat arbitrarily; 
the important thing is that they must be large enough so that the 
filter will forget the Initial values as more data arrive. 

For this linear stochastic control problem, the state-transition 
matrix must be computed carefully to avoid introducing further errors 
into the model equations; therefore, in the study, second-order terms 
are also taken into the computation of Since the state variables 
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of the rainfall-runoff system vary relatively slowly with time, a 
sub-interval of the order of a few minutes is sufficient for the system 
to be considered as time-invariant during this period of time. Thus, 
the interval between two consecutive observations (1 hour) is subdivided 
into only 10 sub-intervals to avoid excessive computational requirements, 
and the value of $ is computed from the approximation formulae (4-20) 
and (4-21) . 

Based on nominal values of the state, the control gains are 
computed first by the Subroutine ADAPT. A sample listing of these 
gains is given in Table 4-1. Later, values of the control gains are 
updated as soon as new values of the estimates are obtained. The 
values of the time-invariant parameters in the model converge relatively 
quickly to their optimal estimates after only 10 iterations, as shown 
in Fig. 4-2. The optimal estimates are: 

K. = 19.99 , ~ = 0.16 , N = 1.18 

1 

Thus, using the same coefficients as those used by Labadie (1968), one 
has K N 

K = 3.77 , A. = ” = 0.16 , N = 1.18 . 

1 Kj ^ h 

These values are not very much different from those obtained by Prasad, 
using numerical integration of the nonlinear equation along with a 
trial and error procedure; Prasad obtained the values 3.79, 0.076 and 
1.27 for A^, A^ and N, respectively. Labadie, using the quasi- 
linearization technique, obtained 4.473, 0.0943 and 1.27, respectively. 
The differences are mainly due to noise terms introduced into the model 
equations to make them more realistic and conform with the nature of the 
rainfall-runoff data. 



CONTROL GAINS 


.539454116E-01 

.534591990E-01 

.529202102E-01 

.523227319E-01 

.516604359E-01 

.509263138E-01 

.501126047E-01 

.492107153E-01 

.482111331E-01 

.471033289E-01 

.458756508E-01 

.445152074E-01 

.430077393E-01 

.413374776E-01 

.394869901E-01 

.374370114E-01 

.351662582E-01 

.326512274E-01 

.298659757E-01 

.267818812E-01 

.233673834E-01 

.195877036E-01 

.154045420E-01 

.107757534E-01 

,565500788E-02 

0 . 


.668759662E-01 

.668278519E-01 

.667745150E-01 

,667153902E-01 

.666498512E~01 

.665772045E-01 

»664966820E-01 

.664074335E-01 

.663085175E-01 

.661988922E-01 

.660774046E-01 

.659427788E-01 

.657936038E-01 

.656283193E-01 

.654452000E-01 

.652423397E-01 

,650176321E-01 

.647687515E-01 

.644931306E-01 

.641879369E-01 

.638500475E-01 

.634760206E-01 

.630620663E-01 

.626040140E-01 

.620972797E-01 

.615384615E-01 


-.408799974E-05 

-.225537746E-05 

.667714257E-07 

.195283176E-05 

.196881815E-05 

.180376887E-05 

.164561338E-05 

.149424528E-05 

.134958663E-05 

.121158868E-05 

.108023272E-05 

.955531026E-06 

.837527969E-06 

.726301129E-06 

.621962609E-06 

.524660405E-06 

.434579901E-06 

.351945452E«06 

.277022066E-06 

.210117154E-06 

.151582362E-06 

.101815444E-06 

,612621616E“07 

.304181887E-07 

.983088953E-08 

0 . 


.819002790E-04 

.411390483E-04 

-.577832919E-05 

-.424874495E-04 

-.390819943E-04 

-.358057368E-04 

-.326663210E-04 

“.296616348E-04 

-.267901305E-04 

-.240508404E-04 

-.214433929E-04 

-.189680315E-04 

-.166256367E-04 

-.144177481E-04 

-.123465907E-04 

-.104151020E-04 

-.862696154E-05 

-.698662243E-05 

-.5499344706-05 

-.417123013E-05 

-.300925850E-05 

-.202132480E-05 

-.121627698E-05 

-.603953703E-06 

-.195219667E-06 

0 . 


-.252710147E-04 

.163041946E-03 

-.138525250E-04 

-.145489795E-03 

-.140989337E-03 

-.129816408E-03 

-.119019061E-03 

-.108598278E-03 

-.985568599E-04 

-.888995269E-04 

-.796330341E-04 

-.707662904E-04 

-.623104887E-04 

-.542792452E-04 

-.466837482E-04 

-.395579165E-04 

-.329085672E-04 

-.267655931E-04 

-.211571493E-04 

-.161148474E-04 

-.116739587E-04 

-.787362217E-05 

-.475705816E-05 

-.237178397E-05 

-.769822890E-06 

0 . 


Table 4-1; 


Sample values of control gains. 
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Values of the estimated surface runoff compared with observed 
runoff are given in Figure 4~3 for control and without control of 
the input disturbance. A sample listing of the results for these 
two cases is given in Table 4-2, in which the peak values of the 
estimated and observed hydrographs are Indicated by small arrows. 

From these results, the following remarks can be made: 

(1) Adaptive control of the inputs is important for the analysis 
of the rainfall-runoff process. Without controlling the rainfall data, 
the system identification results would be bad, and, therefore, the 
estimation of the surface runoff from these noisy inputs would be 
unacceptable . 

(2) Controlling the inputs re<juires additional computer time* 

For the particular data set under study, with the same initial conditions 
mentioned previously, only 6 secs, are needed to add to 12 secs, which 
is the computer time used by the Kalman filtering scheme, including the 
adaptive estimation of the error-covariance matrices R and Q. 

(3) Even with input-control procedure, the approximation of a 
nonlinear system by linearized equations cannot offer good results, 
unless the system under study is not highly nonlinear. 

The effect of adaptive estimation of the model— error covariance 
matrix can be seen from Figure 4-4. In this test case, the Q matrix 
can be seen as 

Q = 






66 


ESTIMATED 


.991807573E-07 
.421776526E-03 
.134001409E-02 
.197257266E-02 
.206403744E-02 
^ .218885424E-02 
.216242101E-02 
.210813662E-02 
.195767699E-02 
.177955574E-02 
.160873416E-02 
.137013063E-02 
.113689880E-02 
.941007460E-03 
.777018873E-03 
.652996927E-03 
.521249447E-03 
.426869681E-03 
.374135431E-03 
*310153269E“03 
^ .240757733E-03 
.184397605E-03 
.133759896E-03 
.100924958E~03 
.659088663E-04 
.298520185E-04 

(I) 


Table 4-2; 


ESTIMATED 


•991807573E-07 

.252270195E-03 

.112551851E-02 

.158049358E-02 

.169617079E-02 

.171930534E-02 

.168333312E-02 

.159430782E-02 

.144564441E-02 

.131368394E-02 

.120310747E-02 

.102278394E-02 

.879052731E-03 

.771580819E-03 

.679555547E-03 

.613792339E-03 

.522779863E-03 

.473724189E-03 

.449590865E-03 

.396322809E-03 

.343625931E-03 

.306895386E-03 

.270461470E-03 

.249638757E-03 

.219309945E-03 

.189475272E-03 

(II) 


MEASURED 


“ 0 . 

.150000000E-03 
.147000000E-02 
.205000000E-02 
.230000000E-02 
*> .240000000E-02 
.238000000E-02 
.225000000E-02 
.200000000E-02 
.178000000E-02 
.160000000E-02 
,128000000E-02 
.103000000E-02 
.850000000E-03 
.700000000E-03 
.600000000E-03 
.450000000E-03 
.380000000E-03 
.360000000E-03 
.280000000E-03 
.200000000E-03 
.150000000E-03 
.lOOOOOOOOE-03 
.800000000E-04 
. 400000000E-04 
- 0 . 


Sample results of estimated surface runoff by 
linear adaptive control approach: 


(I) with control 
(II) without control 
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Without adaptive estimation of the Q matrix at each stage, the filter 
started to diverge at the 20th observation. 

Finally, performances of the adaptive control algorithm with and 
without rectification of the nominal state at each stage are given in 
Figure 4-5 for comparison. A description of subroutine ADAPT and the 
related working subroutines are given in Appendix B. 

Nonlinear Estimation Approach 


Since rainfall data are noisy, the continuity equation is written 
as 

^ = R* - Q + w(t) (4-26) 

where R* denotes the actual noisy input data and w(t) represents the 
input noise which is assumed, as usual, to be a zero-mean white noise 
process with covariance matrix 

E ^ w(t)w(x)'^ j = W 6 (t - t) . (4-27) 

Combining Eqn. (4-26) with the nonlinear storage equation (4-8) , one 
obtains 


Let = Q 


dt^ 


= -(UK - Q) + (1 )„ 

K, K, K„ 



(4-28) 


one then gets the following state equation 
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““ — 


— 





0 

^2 


X -1 ^ 




-X^X3X^X^^ X^ + X^(R* - X^) 



x„ 


0 



3 


+ 

0 

^4 


0 


0 

^5 


0 


0 

L. _ 


*• — 

L J 


(4-29) 

or, in abbreviated notations, 


X = f [x(t),R* ] + g f X(t),w(t) j . 

The observation equation is the same as in the previous section. 
The matrix of partial derivatives is 


(4-30) 


F(t) = 


^2 ^ 

0 

0 

0 


5-1 


X 

-x_x.x_x. 

2 4 5 1 

0 

0 

0 


5-1 


where 


\ = - X2X3X^X^(X3- DX^^ - X^ 




(4-31) 


X5-I 


E3 = - X2X3X^X^ (1 + X3Log(X3)) 


Subroutine ITERA was developed using the iterated extended Kalman 
filter to estimate the state and unknown parameters in the system. 
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Assume that the error-covariance matrices have the following 
values 

R = .001 , W = .01 

and the given initial conditions for X^(0) and P(0) are the same 
as in the previous example. The optimum values of the model para- 
meters converge a little faster than in the case of using the adaptive 
control algorithm; this is shown in Figure 4-6. These optimum estimates 
are 

K, = 19.998 , ” = 0.162 , N = 1.182. 

1 

One could expect that these results are better than those obtained 
previously, since in this case a nonlinear filter has been used and 
therefore model-error has been reduced. This fact is verified by 
a better fit of the estimated outflow with the measured outflow for 
this case as shown on Figure 4-7. A sample listing of the values of 
these estimates is given in Table 4-3. Since the estimates of model 
coefficients converge to stable values, one may conclude that these 
coefficients are constant or can be approximated by time-invariant 
parameters for a particular hydrograph. 

The variation of the estimation-error variances of the state 
with time is shown in Figure 4-8. The effect of adaptive estimation 
of the model-error covariance matrix is also tested in this case. 

Using the same set of initial conditions as used previously, the 
divergence of the filter was less rapid than in the adaptive control 
approach. The result is shown in Figure 4-7. 

A description of subroutine ITERA and related working subroutines 
are also presented in Appendix B. 
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ESTIMATED 


MEASURED 


.119619071E-15 
.150062718E-03 
.146647386E-02 
.205084022E-02 
.221341347E-02 
^ .240452681E-02 
,236577599E-02 
.224040644E-02 
.200300775E-02 
.180874269E-02 
.161181444E-02 
.127554449E-02 
.103091575E-02 
.849778807E-03 
.700038029E-03 
.599990928E-03 
.450001643E-03 
.379999577E-03 
.359997687E-03 
•280013859E-03 
.199952195E-03 
.150323660E-03 
.994354885E-04 
.818362730E-04 
,364464120E-04 
.431316959E-04 


- 0 . 

.150000000E-03 
.147000000E-02 
.205000000E-02 
.230000000E-02 
^ .240000000E-02 
.238000000E-02 
.225000000E-02 
.200000000E-02 
.178000000E-02 
.160000000E-02 
.128000000E-02 
.103000000E-02 
.850000000E-03 
.700000000E-03 
. 600000000E-03 
.450000000E-03 
.380000000E-03 
.360000000E-03 
.280000000E-03 
.200000000E-03 
.150000000E-03 
.lOOOOOOOOE-03 
.800000000E-04 
.400000000E-04 
- 0 . 


Table 4-3: Sample results of estimated surface runoff by nonlinear filter 
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// observation 


Figure 4-8: Variations of the estimation-error variances of the estimates 
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CHAPTER 5 

CONCLUSIONS AND RECOMMENDATIONS 

Summary of the work described in previous chapters 

After reviewing various parameter identification techniques 
used in hydrology in the past, two new approaches for the analysis 
of the rainfall-runoff process are proposed, namely the adaptive 
control of input disturbances and the nonlinear estimation of 
unknown parameters in a noisy environment. Both techniques are 
then applied to the identification of time-invariant parameters 
in the Prasad model of rainfall-runoff. The results obtained 
are encouraging and conform with previous studies made by Prasad 
(1967) and Labadle (1968). 

Advantages of the proposed approaches 

The identification techniques presented in this report 
offer the following advantages: 

1) The formulation of the parameter identification problem 
is natural and gives more insight to researchers about the 
process under study, 

2) Both techniques offer a better and systematic way to 
analyze the rainfall-runoff process. 

3) The identification schemes are sequential and adaptive 
and therefore the proposed approaches can handle large data 

set which may be imbedded in noise with unknown characteristics. 
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4) Both techniques can handle any type of parameters, either 
time -invariant or time dependent . 

5) They can handle multiple-input - multiple-output systems. 

6) The computational requirements are relatively small, since 
the algorithms are simple and the computer times needed to process 
a set of 26 observations are just 14 sec. and 18 sec. for the 
nonlinear estimation approach and the adaptive control approach, 
respectively. 

LIMITATIONS : 

Any approach must have certain operational limitations which stem 
from the formulation of the problem and assumptions used to derive 
necessary equations for further considerations. Therefore, the 
following limitations on the use of the proposed techniques can be 
listed: 

1) They are restricted to jointly Gaussian variables with 
white noise processes or with a certain class of colored noise, i.e., 
noise v(t) which is not white, but can be expressed as 

v(t) = A(t) v(t) + C(t) 

where A(t) is known and C(t) is white. 

2) The system under study must be controllable and observable. 
This problem of stochastic controllability and observability has been 
discussed in detail by Aoki (1967) and Sorenson (1968) . 
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3) When using linearized equations for the adaptive control 
scheme, the initial guess for the value of the state cannot be far 
away from the true value* if filter divergence is to be avoided. 

Suggestions for Future Research 

Possible extensions of the study reported in this report are 
suggested as follows; 

1) Extend this work to the study of distributed-parameter response 
models of hydrologic systems. In this case, the systems will be 
represented by partial differential equations and the estimation 
algorithms developed in Chapter 3 will need only slight modifications 

to include the spatial variations of the states. 

2) Investigate the use of the proposed approaches for short-term 
streamflow prediction in small watersheds, using only measured runoff 
at the mouth of the basins. In this case, the adaptive control scheme 
will be used to estimate the unknown input from measured runoff and 
good short-term streamflow prediction can be obtained by propagating 
the state and the estimation-error covariance matrix forward in time, 

3) Investigate the possibility of combining the identification 
of lumped -parameter models of various environmental zones of a large 
mountaineous watershed in one task, using the proposed identification 
techniques and total measured runoff. In this case, routing models 
for two adjacent zones must be incorporated into the system equations 
and the use of the state-variable approach can give a simple matrix 
differential equation representing the response model of the whole 
watershed under study. Thus, using total measured runoff, one can 
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estimate all unknown parameters in various environmental zones and 
in the routing models at the same time. 

4) Apply these techniques to the analysis of other hydrologic 
processes which can be represented by lumped-parameter response models. 
For this task, no further modifications of the proposed approaches 
are required because they are derived for use in the general case 
of time- invariant or time-dependent, linear or nonlinear lumped- 
parameter hydrologic response models. 
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APPENDIX A 

PROOF OF THE SEPARATION PRINCIPLE FOR 
THE LINEAR OPTIMAL STOCHASTIC CONTROL PROBLEM 



89 


Problem: 


Given the model 


6x = $ , 6x + P _ 6u + w , 

n n,n-l n-1 n,n-l n-1 n-1 


6y = H 6x + V 
n n n n 


(A-l) 


(A-2) 


where 


=0 . E 5 5x^= , 

0 for all n , 

lww'^(= , e\vv"^/=R6 

/n. m) nnm )nm^ nnm 


The various random variables are also assumed to be mutually 
uncorrelated, so 

E I ^ = E j fix^v^ ^ = 0 for all n 

and E < w v^ ( = 0 for all n, m. 

} Tim] 


for all n, m. 


Choose the N control vectors 6u^ (i=0,l, 2, . . , ,N-1) as functions of 

the measurements <Sy^ , 6y^, . . , , 6y^ so that the performance index 


^ r N j 


'(A-3) 


is minimized 



90 


Separation Principle ; 

For the model described by Eqns . (A-1) , (A-2) and (A-3) , the 
optimal stochastic control law is described by 


6u. 


-k-1 “ “ \-k ^N-k,N-k-l 


(A-4) 


where bhe control gain, obtained by solving the deterministic 

control problem, represents the optimal linear estimate of 

the state obtained from the set of measurements 6Y^ k 1’ obtaining 
the estimate, treated as a deterministic quantity. 

Using the definition of expectation, Eqn. (A-3) can be written as 

' 1"n)= J g 


J 


+ <^VnS + 


(A-5) 

where J denotes multiple integrals and 

= d(6x^)d(6xp...d(6xjj)d(6y^)d(6yj^)...d(6y|j_p 

But by hypothesis, the control 6u^ is determined from the measure- 
ments 6Y^, so the and can be integrated out of the first 

term of Eqn. (A-5) 

= fi: + 'Vl®3-l«“j-l>P<«Vl»«V2>''<«^-l>«V2> 

J j=l 



9J 


Consider the last control interval, only the last term of Eqn. (A-6) 
is involved. Using the property of joint probability density function, 
one has 


and 

p ( 6Xjj , ” J P ^ I ®N-1 ’”n-1 ’ ^ ^"n-I I '^^-1 ’ '^^N-1^ 

But the set of measurements determines 6Ujj_^, and 

and 6u^ ^ completely define by Eqn. (A-1) , Thus 

p(6Xjj|6x^_^,w^_l, 6Uj^_l) = 6(6x^- FNjN-I^^-I'^^N-I^ 

(A-9) 


where i5(,) is the Dirac-delta function. 

Since w, is uncorrelated with 6x. and 6y . » one has 
X 1 i 


p^Vil^Vi^^Vi^ = p^-i^ 

Using Eqns. (A-8) , (A-9) and (A-10) , E be written 


^ 1‘^lj “ J ^^-1®N-1'^^-1^ 

. 6(6x^ " '^n,N- 1'^^-1 “ rN,N-l'^^-l“^N-l^ 
• P ^"n- 1> P ^ ’ '^N-l> " ^ 

Integrate Eqn. (A- 11) with respect to 6x^ and then 



9Z 






+ P (A-12) 

only occurs explicitly in (A-12) , so integration will eliminate 

9 ’ control is to be computed as a function of the 6Y , - , so 

N-1 


write 




and choose minimize 


ft 


j ['^'"n- 1^^N-1 rN,N-l\rN,N-l^'^^-l ■'■ ^‘^^-A,N-i\ rN,N-l^^-l 

■p^'VihW'^Vi 

~ '^“n- 1^®N-1 r N,N-l\ Pn,N-1^*“n-1 

^ ' ^H-1 ^ Vi) *n,n-A Tn.n-i^ Vx (A-13) 

The optimal control 6u^ ^ is 

A ^ 

*^^-1 " ~^®N-1 rN,N-l^ Tn^N-I^ ^^N.N-l^^N.N-l"^^-! (A-14) 

^_X “ ^{“^^N-ll '^\-l|' minimum variance estimate of 6x^ 

^ " ^®N-1 ‘‘‘ rN,N-l\ TnjN-I^ (A-15) 

Eqn. (A-14) then becomes 


where 6 
Let 


^'^I-l ■ ^N.N-l ^^-1 


(A-16) 
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Thus, the separation principle is true for the last stage. 

The control (A-16) allows (A-11) to be rewritten as 

E rN,N-l V n,N-1'^^-1 

Let “ ^^- 1 ~ ^^-1 ^^P^^sent the error in the estimate, Eqn. 

(A-17) can then be rewritten in simpler form as 

f[*^^“l^N,N-l\^N,N-l‘^^-l'^ ^^-l^N,N-ArN,N-lVN,N-l'^^-l] 

J \ •* 

Using again the property of joint probability density function, one 
gets 

^ ^^,N-1*’n,N-1^4 Tn.n-i^j 

"u-x*=%-r Vi^p^Vi^ 

•p ( . 6 Yj,_ 2 > d (v^_l , 6 X^_ 2 . 6 Yj,_P ( A-19) 

The first term of the integrand does not depend upon ^ and the error 
is not affected by the control 6 u^ 2 > so 
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i * J 


T 


(A-20) 


where 




Now, consider the last two control intervals 


^{•^ 2 } = «V2V2^“N-2>P(«Vr«\-2>“(«'Sj-l’«V2^ 

The principle of optimality allows one to use expression (A-16) for 
Thus, with Eqn. (A-20), the above equation becomes 


f 


'*'n,N-1*N*N,N-1^'^*N-1‘*' '^“n-2®N-2'^“n-2 


•P<'Vr%-2>‘><Vl’^V2) 

+ rN.N-l\l*N,N-A-l+ Vn-1> <^-22) 

Proceeding in a manner analogous to that employed in deriving 
Eqns. (A-11) - (A-13) , one obtains 


^^-2 “ “ VlVl , N - 2‘^^-2 

where is determined by 

A = / T T p . B ^ T ' 

N-l An-1,N-2 N-1 In- 1,N-2 N-2'^ rN-l,N-2N-l 


(A-23) 


(A-24) 


Vl = 




(A-25) 
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— ^ 

^N-1 " ^-1 ‘ N-2^-1 


The cost associated with 


<S^_2 


is found to be 


^ 1 2 J r^-2^N-l,N-2S-l^N-l,N-2'^^-2^^^-2’'^N-3^^^'^-2’'^^N-3^ 


+ trace ^N,N-A ^ N,N-l\%,N“ A-1 
T ^ 

^N-1,N-2S?-1 -I^N-1 ,N-2\-l\-l ,N-2^N-2 


^N-1 ■'■ ^N-l^N-2 


(A-27) 


The proof of the general result is done inductively. For any k 
one can assume that 


“^Sl-k-l N-k%-k,N-k-l'^^-k“l 

T ^ ~1 

^-k " ^ rN-k,N“k-l^-krN^k,N-k-l ‘‘‘ ®N-k-l^ 


’ ^N-kjN-k-lS-k 


(A-28) 


(A-29) 


^N~k “ "^N-k+ljN-khl-k+l^-k+ljN-k \-k 

I t 

^N-k " h^-k " S«-krN-k,N-k-l^J-k 


(A-30) 


(A-31) 


The cost at stage k is 


T T 



"^^-k^N-k+1 ,N-k^N-k+l^N-k+l jN-k'^^-k^ ^ *^^-k* “^^N-k-l^ 

*^^^\-k’^\-k-l^ 2 Vj+l,N-jS-j+irN-j+l,N~j 
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The proof for the (k+l)st stage is accomplished in the same way 
so that details shall be omitted, Eqns. (A-28) *- (A-32) define the 
optimal stochastic control policy. This completes the proof of the 
Separation Principle. 
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APPENDIX B 
PROGRAM DESCRIPTION 
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SUBROUTINE ADAPT (MX, MEAS, DELT, NSTEPS, AJ, BJ, OX, R, Q, CP, CONTR, Z) 
** This subroutine estimates the state-variables of a linear system 
noisy observations using linear stochastic control approach. The 
optimal estimates are given by a Kalman filter with adaptive estimation 
of the noise covariance matrices ** 

Working parameters : 

MX = number of state-variables, 

MEAS = number of observations, 

DELT = time-period for sub-intervals, 

NSTEPS = number of sub-intervals, 

AJ(.,.) = weighting matrix for the state- variables, 

BJ = weighting factor for the control input, 

0X(.) = set of initial estimates of the state— variables , 

R = observation-error variance, 

Q(.,.) - model-error covariance matrix, 

CP(.,.) = estimation-error covariance matrix, 

CONTR(.) = set of inputs, 

Z(.) = set of observations. 

Subroutines required : STATE, HAD, ADJUSTR, GAIN, STAEST, ERVAR, FPHI, 

PREDICT, CGAIN, GAMMA, UPCON. 

SUBROUTINE STATE (OX, DELT, CTR) 

** This subroutine computes the next value of the state of the 
system ** 

Working parameters : OX, DELT 

CTR = the control input. 

Subroutines required : None . 
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SUBROUTINE HAD (MX, OXP , H) 

** This subroutine computes the observation matrix ** 

Working parameters : MX 

OXP(.) = set of predicted state-variables, 

H(.) = the observation matrix. 

Subroutines required ; None 

SUBROUTINE ADJUSTR (MX, R, WG, DZ, H, XP, DEL) 

** This subroutine adjusts the value of the observation- error variance ** 

Working parameters : MX, R, H 

WG = observation weight, 

DZ = measurement residual, 

XP(.) = set of predicted state-variables, 

DEL = difference between computed and measured residuals. 

Subroutines required : None . 

SUBROUTINE GAIN (MX, PP , H, R, GK) 

** This subroutine computes the Kalman gains ** 

Working parameters : MX, H, R 

PP(,,.) = predicted estimation-error covariance matrix, 

GK(.) = vector of Kalman gains. 

Subroutines required : None . 

SUBROUTINE STAEST (MX, GK, DEL, XP, X) 

** This subroutine gives the best estimates of the state-variables. ** 

Working parameters : MX, GK, DEL, XP 

X(.) = set of best estimates of the state-variables. 


Subroutines required: None. 
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SUBROUTINE ERVAR (MX, PP, GK, H, R, CP) 

** This subroutine updates the estimation-error covariance matrix ** 
Working parameters : MX, PP, GK, H, R 

CP(.,.) = estimation-error covariance matrix. 

Subroutine required : None . 

SUBROUTINE FPHI (MX, OX, CONTR, BELT, NSTEPS, PHI) 

** This subroutine computes the system state-transition matrix, ** 
Working parameters ; MX, OX, CONTR, BELT, NSTEPS 

PHI(.,.) = the system state-transition matrix. 

Subroutines required : STATE . 

SUBROUTINE PREDICT (MX, PHI, GAM, CONGK, CP, Q, X, XP, PP) 

** This subroutine computes predicted values of the state-variables 
and estimation-error covariance matrix ** 

Working parameters : MX, PHI, GAM, CONGK, CP, Q, X, XP 

PP(«j«) = predicted estimation-error covariance matrix. 
Subroutines required : None . 

SUBROUTINE CGAIN (MX, MEAS, AJ, BJ, GAMS, PHIS, CGK, OL) 

** This subroutine computes the control gains ** 

Working parameters : MX, MEAS, AJ, BJ, GAMS, PHIS, CGK 

OL (.,.,.) = set of control cost matrices. 

Subroutines required : None 

SUBROUTINE GAMMA (MX, OX, PHI, GAM) 

** This subroutine is used to compute the model— error covariance 
matrix Q. ** 

Working parameters : MX, OX, PHI 

GAM(.) = control-input matrix. 
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Subroutines required : None. 

SUBROUTINE UPCON (MX, PHI, GAM, OL, AJ, BJ, CONGK, KK) 

** This subroutine is used to update the control law ** 

Working parameters : MX, PHI, GAM, OL, AJ , BJ, KK 

CONGR(.) = set of control gains. 

Subroutines required : None. 

SUBROUTINE ITERA (MX, MEAS, DELT, NSTEPS, X, R, QEL, CP, CONTR, Z) 

** This subroutine estimates the state- variables of a nonlinear system 
through the iterated extended Kalman filter with adaptive estimation of 
the noise covariance matrices ** 

Working parsmieters : MX, MEAS, DELT, NSTEPS, R, CP 

X(.) = state-variable vector 

QEL = variance element of the model-error covariance matrix. 
Subroutines required : ISTATE, IHAD, lADJUSTR, GAIN, STAEST, ERVAR, 

IFPHI, IPREDICT, IGAMMA. 

SUBROUTINE ISTATE (X, DELT, CONTR) 

** This subroutine computes the next value of the state from model 
equations ** 

Working parameters : X, DELT 

CONTR = the control input . 

Subroutines required : None. 

SUBROUTINE IHAD (MX, XP, H, HSM) 

** This subroutine computes the observation matrix and its partial 
derivatives ** 

Working parameters ; MX, H 

XP(.) = set of predicted state-variables, 

HSM = h(X) = observation-vector. 
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Subroutines required: None. 

SUBROUTINE lADJUSTR (MX, R, WG, Z, HSM, H, OXP , SX, X, DEL) 

** This subroutine adjusts the value of the observation-error 
variance ** 

Working parameters : MX, R, WG, Z, HSM, H, OXP, X, DEL SX(.) = set of 

stored predicted-values of the state-variables. 

Subroutines required : None. 

SUBROUTINE IFPHI (MX. X, CONTR, BELT, NSTEPS, PHI) 

** This subroutine computes the system-state- transition matrix ** 
Working parameters : MX, X, CONTR, BELT, NSTEPS, PHI. 

Subroutines required : ISTATE. 

SUBROUTINE IPREDICT (MX, PHI, GAM, CP, QEL, X, NSTEPS, BELT, CONTR, 

XP, PP) 

** This subroutine computes predicted values of the state-variables 
and estimation-error covariance matrix ** 

Working parameters : MX, PHI, GAM, CP, QEL, X, NSTEPS, BELT, CONTR, 

XP, PP. 

Subroutines required : None. 

SUBROUTINE IGAMMA (MX, X, GAM) 

** This subroutine is used to compute the model-error covariance matrix 

Q. ** 

Working parameters : MX, X, GAM. 

Subroutines required : None. 

SUBROUTINE ADJUSTQ (MX, PHI, H, CP, R, SN, DEL, WG, Q) 

** This subroutine adjusts the value of the model— error covariance 
matrix Q ** 
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Working parameters ; MX, PHI, H, CP, R, SN, DEL, WG, Q 
Subroutines required: None. 



